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NUMERICAL  TECHNIQUES  IN  MATHEMATICAL  PROGRAMMING 

t»y 

R.  H.  Bartels,  Q.  H.  Golub,  M.  A.  Saunders 

Abstract 

The  application  of  numerically  stable  matrix  decompositions  to 
minimization  problems  involving  constraints  Is  discussed  and 

shown  to  be  feasible  without  undue  loss  of  efficiency. 

Part  A  describes  computation  and  updating  of  the  product -form 
of  the  UJ  decomposition  of  a  matrix  and  shows  it  can  be  applied  to 
solving  linear  systems  at  least  as  efficiently  as  standard  techniques 
using  the  product-form  of  the  inverse. 

Part  B  discusses  orthogonal! zation  via  Householder  transformations, 
with  applications  to  least  squares  and  quadratic  programming  algorithms 
based  on  the  principal  pivoting  method  of  Cottle  and  Dantzig. 

Part  C  applies  the  singular  value  decomposition  to  the  nonlinear 
least  squares  problem  and  discusses  related  eigenvalue  problems. 
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'introduction 


This  paper  describes  the  application  of  numerically  stable  matrix 
decompositions  to  minimization  problems  involving  linear  constraints. 
Algorithms  for  solving  such  problems  are  fundamentally  techniques  for 
the  solution  of  selected  systans  of  linear  equations,  and  during  the 
last  fifteen  years  there  has  been  a  major  improvement  in  the  understanding 
of  these  and  other  linear  algebraic  problems .  We  show  here  that  methods 
which  have  been  analysed  by  various  workers  and  proven  to  be  numerically 
stable  may  be  employed  in  mathematical  programming  algorithms  without 
undue  loss  of  efficiency. 

Fart  A  describes  means  for  computing  and  updating  the  product-form 
of  the  UJ  decomposition  of  a  matrix.  The  solution  of  systems  of  equations 
by  this  method  is  shown  to  be  stable  and  to  be  at  least  as  efficient 
as  standard  techniques  which  use  the  product -fort:  of  the  inverse. 

In  Fart  B  we  discuss  orthogonalization  via  Householder  transformations. 
Applications  are  given  to  leaB~ 'squares  and  quadratic  programming  algorithms 
based  on  the  principal  pivoting  method  of  Cottle  and  Dantzig  [5  ].  For 
further  applications  of  stable  methods  to  least  squares  and  quadratic 
programming,  reference  should  be  made  to  the  recent  work  of  R.  J.  Hanson  [13] 
and  of  J.  Stoer  [26]  whose  algorithms  are  baaed  on  the  gradient  projection 
method  of  J.  B.  Rosen  [24]. 

In  Part  C  the  application  of  the  singular  value  decomposition  to 
the  nonlinear  least  squares  problem  is  discussed,  along  with  related 
eigenvalue  problems. 


A.  THE  USE  OF  LU  DECOMPOSITION  IN  EXCHANGE  ALGORITHMS 


1.  DU  Decomposition 

If  33  is  an  n  x  n  ,  nonsingular  matrix,  there  exists  a  permutation 
matrix  TT  ,  a  lower-triangular  matrix  L  with  ones  on  the  diagonal,  and 
an  upper-triangular  matrix  U  such  that 

(1.1)  1TB  =  nr  . 

It  is  possible  to  choose  TT  ,  L  ,  and  U  so  that  Ail  elements  of  L 
are  bounded  in  magnitude  by  unity. 

A  frequently-used  algorithm  for  computing  this  decomposition  is 
built  around  GauBSian  elimination  with  row  interchanges.  It  produces 
the  matrices  TT  and  L  in  an  implicit  tom  as  Bhown: 

For  k  =  1,2,  ...,n-l  in  order  carry  out  the  following 
two  steps: 

4 

(1.2)  Find  an  element  in  the  k-th  column  of  B  ,  on  or  below  the 
diagonal,  which  has  maximal  magnitude-  Interchange  the 
k-th  row  with  the  row  of  the  element  found. 

(1.3)  Add  an  appropriate  multiple  of  the  resulting  k-th  row  to  each 
row  below  the  k-th  in  order  to  create  zeros  below  the 

in  the  k-th  column. 
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Each  execration  of  the  first  step  (1.2),  in  matrix  notation,  amounts 
to  the  premultiplication  of  B  by  a  suitable  permutation  matrix  TT.  . 
The  following  step  (1.3)  may  be  regarded  as  the  premultiplication  of  B 
by  a  matrix  Tk  of  the  form 

k 


(1.4)  k 


where  Ig,  .  I  <  1  for  each  i  *  kt-1, ...,n  . 
x,  K'  — 

By  repeating  the  two  steps  n-1  times,  B  is  transformed  into  U  . 

And  at  the  same  time  the  matrix  (L  "Ht)  is  collected  In  product  form 

(1.5)  r„-lVl--rl"L  • 

x  o 

This  algorithm  requires  nr/ 5  +  0(n  )  multiplication/aivision  operations 
and  again  this  many  addition/ subtraction  operations.  Both  U  and  all 
of  the  gA  j  can  be  stored  in  the  space  which  was  originally  occupied  by  B  . 
An  additional  n  locations  are  required  for  the  essential  information  contained 
in  the  Tf^  . 
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Many  algorithms  require  the  solving  of  a  sequence  of  linear  equations 
)  B^x  =  v^ 


for  which  each  Bv;  differs  front  its  predecessor  in  only  one  column. 
Examples  of  such  algorithms  are:  the  simplex  method,  Stiefel's  exchange 
method  for  finding  a  Chebyshev  solution  to  an  overdetermined  linear 
equation  system,  and  adjacent  -  path  methods  for  solvin':  the  complementary- 
pivot  programming  problem. 

Given  that  has  a  decomposition  of  the  form 

(2.2)  , 

where  t/0^  is  upper-triangular,  and  given  that  has  been 

stored  ae  a  product 


l(°)’1  .  r(0)  -(0)  r(o)ff(o) 

L  *  rn-lVl  ”•  fl  *1  » 


the  initial  system  of  the  sequence  is  readily  solved:  Set 

(2.1.)  , 


and  then  hack-solve  the  triangular  system 


„«»,  =  , 


$ 

fc 


3 .  Updating  the  LU  Decomposition 


Let  the  column  r^  of  be  replaced  by  the  column  vector  . 

So  long  as  we  revise  the  ordering  of  the  unkrowns  accordingly,  we  may 
insert  into  the  last  column  position,  shifting  columns  r  +1 

through  n  of  one  position  to  the  left  to  make  room.  We  will 

fl) 

call  the  result  s  ~  ,  and  wc  can  easily  check  that  it  has  the 
decomposition 

(3.1)  =  L^V1)  , 


/  V\ 

where  H'  1  is  c  matrix  which  is  upper -Hessenberg  in  its  last  n  -  tq +  1 
columns  and  upper-triangular  in  its  first  r^-1  columns.  That  is, 
has  the  form 


The  first  r  -1  columns  of  H 
o 


(1) 


(0) 


are  identical  with  those  of  U 

(0) 


The  next  n-r  are  identical  with  the  last  n-r  columns  of  U 
o  o  > 

And  the  last  column  of  is  the  vector  i/0^  . 


H 


(1) 


can  be  reduced  to  upper-triangular  form  by  Gaussian  elimination 


with  row  interchanges.  Here,  however,  we  need  only  concern  ourselves 
with  the  interchanges  of  pairs  of  adjacent  rows.  Thus 


is  gotten 


from  H'  fay  applying  a  sequence  of  aiaple  trana format iona: 


(1)  (1)  (1)  (1)  (1)  (1) 

u  rn-l”n-l  *•*  rt  “r  H  * 

o  o 


where  each  has  the  fora 


(3A) 


and  each  TT^  ia  either  the  identity  matrix  or  the  identity  with  the  i-th 
and  i+l-st  rows  exchanged,  the  choice  being  made  so  that  |g^|  <  l  . 

The  essential  information  in  all  of  these  transformations  can  be 
stored  in  n-rQ  locations  pins  an  additional  n-re  bits  (to  indicate 
the  interchanges) .  If  we  let 


=  r(1)tr(1) 
‘n-rn-l 


r  r 
o  o 


then  we  have  achieved  the  decomposition 

(3-6)  3^  --  l/1V'±'1 


The  transition  from  E'  •  to  b'  '  for  any  i  is  to  be  made 
exactly  as  was  the  transition  from  to  3'^  .  Any  system  of 

linear  equations  involving  the  matrix  for  any  i  is  to  be  solved 

according  to  the  steps  given  in  (3,4)  and  (2.5). 


I 
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4.  Round-off  Considerations 


For  most  standard  computing  machines  the  errors  in  the  basic 
arithmetic  operations  can  be  expressed  as  follows: 

fl(a  +  b)  =  afl  +  e1)  +  b(l  +  e2) 

(4.1)  fl (a  x  b)  =  ab(l  +  e^) 

«(®A)  »  (a/b)(l  +  , 

where  |e.J  <  .  Here  p  stands  for  the  base  of  the  number  system 

in  which  machine  arithmetic  is  carried  out  and  t  is  the  number  of 
significant  figures  which  the  machine  retains  after  each  operation.  The 
notation  ff (a  "op"  b)  stands  for  the  result  of  the  operation  "op" 
upon  the  two,  normal -precision  floating-point  numbers  a  and  b  when 
standard  floating-point  arithmetic  is  used. 

The  choice  of  an  UJ  decomposition  for  each  and  the  particular 

way  in  which  this  decomposition  is  updated  were  motivated  by  the  desire 
to  find  a  way  of  solving  a  sequence  of  linear  equations  (2.1)  which  would 
retain  a  maximum  of  information  from  one  stage  to  the  next  in  the  sequence 
and  which  would  be  as  little  affected  by  round-off  errors  as  possible. 

Under  the  assumption  that  machine  arithmetic  behaves  as  given  in  (4.1), 
the  processes  described  in  Sections  2  and  5  are  little  affected  by 
round-off  errors.  The  efficiency  of  the  processes  will  vary  from  algorithm 
to  algorithm,  but  we  will  argue  in  a  subsequent  section  that  the  processes 
should  cost  roughly  as  much  as  those  based  upon  product-form  inverses 
of  the  . 
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We  will  now  consider  the  round-off  properties  of  the  basic  steps 


described  in  Sections  2  and  3. 

The  computed  solution  to  the  triangular  system  of  linear  equations 


(U.l) 


=  y 


car,  be  shown,  owing  to  round-off  errors,  to  satisfy  a  perturbed  system 
(lt.3)  (U^  +  6U^)x  =  y  . 


It  is  shown  in  Forsythe  and  Moler  [ 9  ]  that 


(t.«0 


!loa(i)|! 

!!u(i)i! 


<  (1.0l)31"t 


9 


where  l|...j|  denotes  the  infinity  norm  of  a  matrix,  and  thus  round-off 
errors  in  the  back- solution  of  a  triangular  system  of  linear  equations 
may  he  regarded  as  equivalent  to  relatively  small  perturbations  in  the 
original  system. 

Similarly,  the  computed  L  and  U  obtained  by  Gaussian  elimination 
with  row  interchanges  from  an  upper -Hessenberg  matrix  H  satisfy  the 
perturbed  equation 


(14-5)  H  +  6H  =  HJ  , 


where  Forsythe  and  Moler  show  that 


(*.6) 


M  < 

IlHlI  " 


2  Dl-t 
n  pp 


I 


and  Wilkinson  [28]  establishes  that  p  <  n  .  Thus,  the  computational 
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process  indicated  in  (3*3)  can  be  regarded  as  Introducing  only  relatively 
small  perturbations  in  each  of  the  . 

Similar  results  hold  for  the  initial  1U  decomposition  (2.2)  with 
a  different  bound  for  p  .  The  reader  is  referred  again  to  Forsythe 
and  Moler. 


The  most  frequent  computational  step  in  the  processes  which  we  have 
described  is  the  application  of  one  Gaussian  elimination  step  r  to  a 
column  vector  v  s 


i 

(4.7)  w  =  Tv  = 

i 


i  i 


The  computed  vector  w  satisfies 

(4.8)  w^  =  vk  for  k/j 

Wj  -  ff(fl(gvi)  +  Vj) 

=  gvjd  +  c?)(l  +  +  Vj(l  +  Eg) 

=  ey±  +  Vj  +  gvjL(e1  +  Ej  +  s^j)  +  Vj62 
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Thus  we  may  regard  the  computed  vector  w  as  the  exact  result  of  a 
perturbed  transformation 

(4.9)  v  =  (r  +  5r)v  , 

where 

i 

(4.10)  5r  = 

j 


i  a 


an  d 

(4.11)  a  =  g(a1  +  e3  + 

1  =  Eg  • 

Therefore  we  have 

lhJS)  M  <  e?  *  y?!  *  M 

Ill’ll  1+l«l 

where  the  right-hand  side  is  bounded,  since  |g|  <  1  ,  according  to 
(4.13)  —  <  P1‘tt3  +  01_t]  <  3.Ol01_t  (say). 

Urn  ~ 

Hence,  the  computations  which  we  perform  using  transformations  (4.7)  also 
introduce  relatively  small  perturbations  into  the  quantities  which  we  manipulate. 


It  is  precisely  with  regard  to  such  transformations  that  ve  feel 
our  method  of  computation,  based  upon  UJ  decompositions,  is  superior 
to  methods  based  upon  the  inverses  of  the  matrices  .  Such  methods 

use  transformations  of  the  form 


alternatively,  in  product -form  methods,  they  are  applied  to  the  vector 
to  produce  the  solution  to  system  (2.1).  As  such,  they  involve 
successive  computations  of  the  form  (h.7) .  Each  such  computation  may  be 
regarded  as  satisfying  (^.9).  But,  since  the  tjj  may  be  unrestricted  in 
magnitude,  no  bound  such  as  (t.lj)  can  be  fixed. 


5*  Efficiency  Considerations 


Afl  we  nave  already  pointed  out,  it  requires 
(5.1)  n3/3  +  0(n2) 

multiplication-type  operations  to  produce  an  initial  IAJ  decompoBition  (2.2). 
To  produce  the  product-form  inverse  of  an  n  x  n  matrix,  on  the  other 
hand,  requires 

(5-2)  n3/2  +  0(n2) 

operations. 

The  solution  for  any  system  (2.1)  must  be  found  according  to  the  LCJ- 
decoraposition  method  by  computing 

(5-3)  y  *  L(i)  1v(i) 


followed  by  solving 

(5. 4)  U^x  =  y  . 


The  application  of  L 


(0) 


-1 


to  v 


(i) 


(5-5) 


2 


in  (5-5)  will  require 


operations.  The  application  of  the  remaining  transformations  in  1^ 
will  require  at  most 

(5.6)  i(n-l) 

operations.  Solving  (5.4)  costs 


12 


(5-7) 


n(n+l) 

2 

operations.  Hence,  the  cost  of  (5*3)  and  (5.4)  together  ifi  not  greater  than 


(5.8) 


n  +  i(n-l) 


£  i 

operations,  and  a  reasonable  expected  figure  would  be  n  +  ±  (n-1)  . 

On  the  other  Land,  computing  the  solution  to  (2.1)  using  the  usual 
product  form  of  B'  '  requires  the  application  of  n+i  transformations 
of  type  (4. it)  to  v^  at  a  cost  of 

,  ,  o 

<  5  •  9 )  n  +  in 


operations. 

If  a  vector  replaces  column  r.  in  B^  ,  then  the 

(i)-1 

updating  of  B'  requires  that  the  vector 


(5.10) 


z  =  B^ 


be  computed.  This  will  cost  n^+in  operations,  as  shown  in  (5.9).  Then 
a  transformation  of  form  (4. it)  must  be  produced  from  z  ,  and  this  will 
bring  the  total  updating  cost  to 


(5*U) 


n  +  (i  +  l)n  . 


The  corresponding  cost  for  updating  the  LU  decanposition  will  be  not  more 
than 


(5-12) 


+  i(n-l) 


operations  to  find  i/^  ,  followed  by  at  most 
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(5-15) 


Efctil 

2 

operations  to  reduce  i+1^  to  and  generate  the  trane format ions 

of  type  (3 .fc)  which  effect  this  reduction*  This  gives  a  total  of  at  most 

(5.11)  n2  +  i(n-l) 

2  1 

operations,  with  an  expected  figure  closer  to  n  + 1  (n-l)  . 

Hence,  in  every  ease  the  figures  for  the  UJ  decomposition:  (5.1k), 
(5.8),  and  (5-1)  are  smaller  than  the  corresponding  figures  (5. 11),  (5-9), 
and  (5*2)  for  the  product-form  inverse  method. 


6.  Storage  Considerations 

All  computational  step*  for  the  LU-dec opposition  method  may  be 
organized  according  to  the  columns  of  the  matrices  .  For  large 

systems  of  data  this  permits  a  two-level  memory  to  be  used,  with  the 
high-speed  memory  reserved  for  those  columns  being  actively  processed. 

The  organization  of  Gaussian  elimination  by  columns  is  well-known, 
and  it  is  clear  how  the  processes  (5 .3)  may  be  similarly  arranged. 
Finally,  the  upper-triangular  systems  (5.^)  can  be  solved  columnwise 
as  indicated  below  in  the  I*  x  4  case: 


U11  u12  u13  U1M  X1 


0  ^22  ^3  x2 

0  0  U33  x3 


0  0  0  u. 


Bring  the  y  vector  and  the  last  column  of  U  into  high-speed 
mosory.  8et  =  y^/u^  .  Set  «  yi  "uiUxi+  for  1  “  3>2,1 
This  leaves  us  with  the  following  3x3  system: 

I  “u  ”12  “15I  *l\  yl' 

)  0  V  ”2}  X2  -  y2 


Ve  process  it  as  suggested  in  the  U  x  1  case,  using  now  the  third 
column  of  U  to  produce  .  Repeat  as  often  as  necessary. 
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fil 

In  the  event  that  the  matrices  3  are  sparse  as  well  as  large, 
we  wish  to  organize  computations  additionally  in  such  a  way  that  this 
sparseness  is  preserved  as  much  as  possible  in  the  decompositions. 

For  the  initial  decomposition  (2.2),  for  example,  we  would  wish  to 
order  the  columns  of  in  Such  a  way  that  the  production  of  I/0' 

and  i/0^  introduce  as  few  new  nonzero  elements  as  possible.  And  at 
subsequent  stages,  if  there  is  a  choice  in  the  vector  which  is 

to  be  introduced  as  a  new  column  into  the  matrix  to  produce 

it  may  be  desirable  to  make  this  choice  to  some  extent  on  sparseness 
considerations. 

It  is  not  generally  practical  to  demand  a  minimum  growth  of  nonzero 
elements  over  the  entire  process  of  computing  the  initial  decomposition. 
However,  one  can  easily  demand  that,  having  processed  the  first  k-1 
columns  according  to  (1.2)  and  (1-3),  the  next  column  be  chosen  from  those 
remaining  in  such  a  way  as  to  minimize  the  number  of  nonzero  elements 
generated  in  the  next  execution  of  steps  (1.2)  and  (1.3) •  See,  for 
example,  Tewarson  [27]  Choice  of  the  next  column  may  also  be  made 
according  to  various  schemes  of  "merit” ;  e.g.,  see  Dantzig  et  al.  [6]« 

The  Introduction  of  new  nonzero  elements  during  the  process  of 
updating  the  i-th  decomposition  to  the  i+l-st  depends  upon 

(6.3)  the  nonzero  elements  in  over  those  in  , 

and 

(6. 4)  the  number  of  the  column  to  be  removed  from  . 
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once 


IJo  freedom  is  possible  in  the  reduction  of  H'*+^  to 

has  been  chosen  and  the  corresponding  r^  has  been  determined. 

The  growth  (£.J)  can  be  determined  according  to  the  techniques 
outlined  In  Tewarson's  paper*  at  a  cost  for  each  value  of  i  *  however* 
which  is  probably  unacceptable.  The  more  important  consideration  is  (f.h), 
The  larger  the  value  of  r^  *  the  fewer  elimination  stepB  must  be  carried 
out  on  and  the  less  chance  there  is  for  nonzero  elements  to  be 

generated.  Again,  howevsr,  the  determination  of  the  value  of  r^ 
corresponding  to  each  possible  choice  of  may  prove  for  most 

algorithms  to  be  unreasonably  expensive. 


Accuracy  Considerations 


During  the  execution  cf  an  exchange  algorithm  it  sometimes  becomes 

necessary  to  ensure  the  highest  possible  accuracy  for  a  solution  to  one 

of  the  systems  (2.1) .  High  accuracy  is  generally  required  of  the  last 

solution  in  the  sequence,  and  it  may  be  required  at  other  points  in  the 

sequence  when  components  of  the  solution,  or  numbers  computed  from  them, 

approach  critical  values.  For  example,  in  the  simplex  method  inner 

products  are  taken  with  the  •  oetor  of  simplex  multipliers,  obtained  by 

(O 

solving  a  system  involving  B  '  ,  and  each  of  the  non-bacic  vectors. 

The  computed  values  are  then  subtracted  from  appropriate  components  of 
the  cost  vector,  and  the  results  are  compared  to  zero.  Those  which  are 
of  one  sign  have  importance  in  determining  how  the  matrix  is 

to  be  obtained  from  .  The  value  zero,  of  course,  is  critical. 

The  easiest  way  of  ensuring  that  the  computed  solution  tc  a  system 

(7-1)  Bx  =  v 

has  high  accuracy  is  by  employing  the  technique  of  iterative  refinement 
[9  ,  Chapter  13].  According  to  this  technique,  if  is  any  sufficiently 

good  approximation  to  the  solution  of  (7.1)  (for  example,  a  solution 
produced  directly  via  the  LQ- decomposition  of  B  )  then  improvements  may 
be  made  by  computing 

(7.2)  =  v  -  Bx'^  , 
solving 

(7.3)  Bz^  =  r^  , 
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and  setting 

(7 .k)  x(J+l)  =  x(d)  +  z(J) 

for  j  »  0,1,2,...  until  Ijz^lj  is  sufficiently  small.  The  inner 
products  necessary  to  form  the  residuals  (7*2)  must  be  computed  in 
double-precision  arithmetic.  If  this  rule  is  observed,  however,  and  if 
the  condition  of  the  system,  measured  as 

(7*5)  cond(B)  *  |!  B  jj  jlB'1!!  , 

"t-1 

is  not  close  to  0  ,  the  refinement  process  can  be  counted  on  to 

terminate  in  a  few  iterations.  The  final  vector  x^'  will  then  be 
as  accurate  a  solution  to  (7.1)  as  the  significance  of  the  data  in  B 
and  v  warrant . 

Step  (7.5)  is  most  economically  carried  out,  of  course,  via  the 
same  LU-decomposition  which  was  used  to  produce  .  If  this  is 

done,  each  repetition  of  steps  (7.2)  through  (7.U)  will  cost  only 
0(n  )  operations.  The  alternative  approach  of  producing  a  highly 
accurate  solution  to  (7.1)  by  solving  the  system  entirely  in  double¬ 
precision  arithmetic  is  generally  more  expensive  than  iterative 
refinement  by  a  factor  of  n  . 
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B.  THE  QR  DECOMPOSITION  AND  QUADRATIC  PROGRAMMING 


8.  Householder  Trlangularization 


Householder  transformations  have  been  widely  discussed  in  the 
literature.  In  this  section  we  are  concerned  with  their  use  in  reducing 
a  matrix  A  to  upper -triangular  form,  and  in  particular  we  wiah  to  show 

how  to  update  the  decomposition  of  A  when  its  columns  are  changed  one 

\ 

by  one.  This  will  open  the  way  to  the  implementation  of  efficient  and 
stable  algorithms  for  solving  problems  involving  linear  constraints. 

Householder  transformations  are  symmetric  orthogonal  matrices  of 
the  form  P^  =  1  -  where  is  a  vector  and  0^  *=  2/  (u^_u^)  . 

Their  utility  in  this  context  is  due  to  the  fact  that  for  any  non-zero 
vector  a  it  is  possible  to  choose  in  such  a  way  that  the 
transformed  vector  P^a  is  zero  except  for  its  first  element. 

Householder  [ 15]  used  this  property  to  construct  a  sequence  of  transformations 
to  reduce  a  matrix  to  upper-triangular  form.  In  [29],  Wilkinson  describes 
the  process  and  his  error  analysis  shows  it  to  be  very  stable. 

Thus  if  A  =  (a^, ...,an)  is  an  m x n  matrix  of  rank  r  ,  then 
at  the  k-th  stage  of  the  triangularization  (k  <  r)  we  have 


Pk-1  Pk-2  *  ’  *  P0A  = 


R,  S, 

k  k 


0  T,. 


where  R^  is  an  upper-triangular  matrix  of  order  k  .  The  next  step 
is  to  compute  A^1^  =  A'^  where  P^  is  chosen  to  reduce  the  first 


20 


r<l»;  ft 


column  of  to  zero  except  for  the  first  component .  This  component 
becomes  the  last  diagonal  element  of  and  since  its  modulus  is 

equal  to  the  Euclidean  length  of  the  first  column  of  it  should  in 
general  be  maximized  by  a  suitable  interchange  of  the  columns  of 

\  .  After  r  steps.  T ,  will  be  effectively  zero  (the  length 

/ 

of  each  of  its  columns  will  be  smaller  than  some  tolerance)  and  the 
process  stops. 

Hence  we  conclude  that  if  rank(A)  -  r  then  for  some  permutation 
matrix  h  the  Householder  decomposition  (or  "QR  decomposition")  of  A  is 


r  n-r 


where  Q  =  Pr_^  Pr_2  . . .  PQ  is  an  m  x  m  orthogonal  matrix  and  R  is 

upper-triangular  and  non-singular. 

We  are  now  concerned  with  the  manner  in  which  Q  should  be  stored 

and  the  means  by  which  Q  ,  R  ,  S  may  be  updated  if  the  columns  of  A 

are  changed.  We  will  suppose  that  a  column  is  deleted  from  A  and 

that  a  column  a  is  added.  It  will  be  clear  what  is  to  be  done  if  only 
q  ^ 

one  or  the  other  takes  place. 

Compact  Method: 

Since  the  Householder  transformations  are  defined  by  the  vectors 
the  usual  method  is  to  store  the  u^'s  in  the  area  beneath  R  ,  with 
a  few  extra  words  of  memory  being  used  to  store  the  &k's  and  the  diagonal 
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elements  of  R  .  The  product  Qz  for  some  vector  z  is  then  easily 

computed  in  the  form  P  ^  Pr_2  ...  PQ  z  where,  for  example, 

?0Z  =  (I  -  g0u0uQ)z  =  z-g0(u0z)uQ.  The  updating  is  best  accomplished 

as  follows.  The  first  p-1  columns  of  the  new  R  are  the  same  as 

before;  the  other  columns  p  through  n  are  simply  overwritten  by 

columns  a  a  ,a  and  transformed  by  the  product  P  ,  P  -  ...  P. 

p+j.  n  q  p-i.  p-2  0 

J;  then  T  ,  is  triangularized  as  usual. 

Tp-1  / 

This  method  allows  Q  to  be  kept  in  product  form  alwayB,  and  there  is  no 
accumulation  of  errors.  Of  course,  if  p  =  1  the  complete  decomposition 
must  be  re-done  and  since  with  m  >  n  the  work  is  roughly  proportional 
to  {m-n/3)n  this  can  mean  a  lot  of  work.  But  if  p  =  a/2  on  the 
average,  then  only  about  l/8  of  the  original  work  must  be  repeated 
each  updating. 

Explicit  Method: 

The  method  just  given  is  probably  best  when  m  »  n  .  Otherwise 
vc  propose  that  Q  should  be  stored  explicitly  and  that  the  updating 
be  performed  as  follows : 

(l)  The  initial  3  can  be  computed  by  transforming  the  identity 
matrix  thus: 

Pr-lPr-2  "•  V"rlI„>  * 


to  obtain  a  new 
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(2)  If  a  is  added  to  A  then  compute  s  =  Qa  and  add  it 
*  Q  Q. 


to  the  end  of 


(?) 


(5)  Delete  a^  where  applicable  (p  <  r)  .  This  normally  means 
just  updating  the  permutation  vector  used  to  describe  IT  . 

(1*)  The  initial  situation 


Q  ATT 


P 


has  thus  been  changed  to 


QAff  = 


where  the  areas 


©,©,©,© 


are  the  same  as  before. 


23 


This  is  analogous  to  the  Hessenberg  form  encountered  in 

updating  HJ  decompositions.  We  now  employ  a  sequence  of 

(r-p)  plane  rotations,  as  used  by  Givens  and  analysed 

by  Wilkinson  [  30],  to  reduce  the  subdiagonal  of  area 

to  zero.  This  changes  areas  and  >  and  the 

•  corresponding  rows  of  Q  must  also  be  transformed.  Since 

the  plane  rotations  are  elementary  orthogonal  transformations, 

* 

the  latter  step  produces  a  new  matrix  Q  which  is  also 

orthogonal,  and  the  work  necessary  is  approximately  proportional 

2 

to  2mn  +  n  . 

(5)  Finally,  a  single  Householder  transformation  Py  is  applied 
_  * 

to  produce  Q  =  ,  where  this  transformation  is  the  one 

which  reduces  area  (?)  to  zeros  except  for  the  first 
element.  The  work  involved  is  proportional  to  2(m-n)m  . 

Thus  the  transformation  Q  reduces  ATT  to  a  new  upper-triangular 
form,  and  the  original  transformations  PQ, • • •, ^  ,  the  plane  rotations, 
and  the  final  Householder  transformation  may  all  be  discarded  since  the 
required  in  format  ion  is  all  stored  in  Q  .  The  total  work  involved  is 

o  p  g 

roughly  proportional  to  (2mn  +  n  )  +  2(m  -  n)m  *  2m"  +  n  and  the  stability 
of  the  orthogonal  transformations  is  such  that  accumulation  of  rounding 
errors  during  repeated  applications  of  the  updating  process  should  be 
very  slight. 


2k 


9-  Project Iona 


i- 

t 

■i. 

i 


£ 


In  optimization  problems  involving  linear  constraints  It  is  often 
necessary  to  compute  the  projections  of  same  vector  either  into  or 
orthogonal  to  the  space  defined  by  a  subset  of  the  constraints  (usually 
the  current  "basis") .  In  this  section  we  shew  how  Householder 
transformations  may  be  used  to  compute  such  projections.  As  we  have 
shown,  it  is  possible  to  update  the  Householder  decomposition  of  a 
matrix  when  the  number  of  columns  in  the  matrix  is  changed,  and  thus 
ve  will  have  an  efficient  and  Btable  means  of  orthogonalizing  vectors 
with  respect  to  basis  sets  whose  component  vectors  are  changing  one  by 
one. 

Let  the  basis  set  of  vectors  a^,  a^,  *  *  *  >  an  form  the  columns  of 
an  m  x  n  matrix  A  ,  and  let  Sy  be  the  sub-space  spanned  by  {a^  . 
We  shall  assume  that  the  first  r  vectors  are  linearly  independent 
and  that  rank(A)  =  r  .  In  general,  m  >  n  >  r  ,  although  the  following 
is  true  even  if  m  <  n  . 

Given  an  arbitrary  vector  z  we  wish  to  compute  the  projections 


u  »  h  ,  v  =  ( I  -  P)  z 


for  seme  projection  matrix  P  ,  such  that 

(a) 

Z  *  U  +  V 

(b) 

T 

u  v  *  0 

(c) 

ueSr  (i.e.,  gx  such  that  Ax  * 

u) 

(d) 

v  ie  orthogonal  to  Sy  (i.e.. 

T 

A  v 
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“f  + 

One  method  is  to  write  P  as  AA  where  A  is  the  n  x  m  generalized 
inverse  of  A  ,  and  in  C  7  3  Fletcher  shows  how  A+  may  be  updated 
upon  changes  of  basis.  In  contrast,  the  method  based  on  Householder 
transformations  does  not  deal  with  A+  explicitly  but  instead  keeps 
AA+  in  factorized  form  and  simply  updates  the  orthogonal  matrix  required 
to  produce  this  form.  Apart  from  being  more  stable  and  just  as  efficient, 
the  method  has  the  added  advantage  that  there  are  alwiys  two  orthonormal 
sets  of  vectors  available,  one  spanning  Sy  and  the  other  spanning  its 
complement . 

As  already  shown,  we  can  construct  an  m  x  n  orthogonal  matrix  Q 
such  that 


r  n-r 


where  R  is  an  rxr  upper-trlar jular  matrix.  Let 


(9*1)  w  =  Qz 

and  define 


3  r 
}  m-r 


(9*2) 


Then  it  is  easily  verified  that  u,v  are  the  required  projections  of  z  , 
which  is  to  say  they  satisfy  the  above  four  properties.  Also,  the  x  in 
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(c)  is  readily  shown  to  he 


fjf  V 


x  = 


In  effect,  we  are  representing  the  projection  matrices  in  the  form 


(9*5) 


„  I 

P  -  ft1!  r  J(Ir  0)Q 

0  ' 


ana 


(9-4) 


I  -  p  =  ft 


»  K..r y- 


m-r  J 


and  we  are  computing  u  =  Pz  ,  v  =  (i-P)z  by  means  of  (9.1),  (9-2). 
The  first  r  columns  of  ft  span  Sr  and  the  remaining  m-r  span 
its  complement .  Since  ft  and  R  may  be  updated  accurately  and 
efficiently  if  they  are  computed  using  Householder  transformations,  we 
have  as  claimed  the  means  of  orthogonalizing  vectors  with  respect  to 
varying  bases. 

As  an  example  of  the  use  of  the  projection  (9*4),  consider  the 

T  T 

problem  of  finding  the  stationary  values  of  x  Ax  subject  to  x  x  =  1 
and  CTx  =  0  ,  where  A  is  a  real  symmetric  matrix  of  order  n  and  C 
is  an  nxp  matrix  of  rank  r  ,  with  r  <  p  <  n  .  It  is  shown  in  [12] 
that  if  the  usual  Householder  decomposition  of  C  is 


r  n-r 
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then  the  problem  li  equivalent  to  that  of  finding  the  eigenvalues  and 

* 

eigenvectors  of  the  matrix  PA  >  where 

/°  0  *\ 

p  -  I  -  P  -  Q  [  R 

V°  W 

is  the  projection  matrix  in  (9.1*) .  It  can  then  be  shown  that  if 


where  ia  r  x  r  ,  then  the  eigenvaluee  of  PA  are  the  same  as 

those  of  and  so  the  eigensystem  has  effectively  been  deflated 

by  the  number  of  independent  linear  constraints.  Similar  transformations 

T 

can  be  applied  if  the  quadratic  constraint  is  x  Bx  =  1  for  Borne  real 
positive  definite  matrix  B  . 
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10.  Orthogonalizatlon  wit;-,  hesceot  to  jcsitiy  '.efirite  Forms 

Flutcner  also  sh  owe  ; ;i  !  .  i  how  to  update  prove  -•!. Lon  matrices  w 
Lv  required  to  orthogonalise  with  respect.  to  a  ,-i/on  posit  I  .e 
definite  matrix  '■  .  now  shew  ;n >w  to  compute  such  prc.icct  i ons  us 
r:ousi*r- o.bier  \  rn.ns  i  ?:\v. ,  a.*ri  dc^ci1*  t.ho  corrjricr. *rs  nia:.io  in  “the  lad 


c?0  *  •.  J ' !  C*  •  *  •  *n 


r>r'.  o 


«wm1m  ;,*a>  »/h*  SJ^pi.’irrti  fivViT  • 

an  arbitrary  z  it  is  required  to  find  n  -=  Pz  , 

T  *  P)  -  for  some  ?'  ,  rue:  h  t  ho  t 


.  S'. , 


i  +  V 


uADv  ■  o 

(c)  lx  suon  tat. s,  Ax  =  u 

(d)  (r-A)1'/  <  o  . 

for  simplicity  we  will  assume  that.  ra-duA)  n  .  Then,  rather  than 
computing  F  explicitly  as  Fletcner  does  according  to 

P  -  A'ArDA)_1  Arn  , 

we  obtain  the  Cholesky  decomposition  of  D  thus: 

■p 

n  =  ll1 


where  L  is  lower -triangular  and  non-singular  if  D  is  positive 

T 

definite.  We  then  compute  B  =  L  A  and  obtain  the  decomposition 


Defining 


v 


and 


u 


}n 

}  m-n 


it  is  easily  verified  that  u, v  are  the  required  projections,  and 

again  the  x  in  (c)  is  given  by  x  =  R-^  .  Since  changing  a  column 

T 

of  A  is  equivalent  to  changing  the  column  L  of  B  ,  the 
matrices  Q  and  R  may  be  updated  almost  as  simply  as  before. 


We  first  consider  minimisation  of  quadratic  forms  subject  to 
linear  equality  constraints.  The  solution  is  given  by  a  single  system 
cf  equations  and  the  algorithm  we  describe  for  solving  tnis  system  will 
serve  as  a  basic  tool  for  solving  problems  with  inequality  constraints, 
it  will  also  provide  an  example  of  how  solutions  to  even  strongly 
ill-conditioned  problems  may  be  obtained  accurately  if  orthogonal inert  ion 
techniques  are  used. 

Let  AjG  be  given  matrices  of  orders  m  x  n  ,  pxn  respectively 
and  let  b,h  be  given  vectors  of  consistent  dimension.  The  least 
squares  problem  to  be  considered  here  is 

Problem  LS:  min  !;b  -  Axj!2 

subject  to  Gx  -  h  . 

Similarly,  let  D  be  a  given  positive  semi -definite  matrix  and  c 
a  given  n-dimensional  vector.  The  quadratic  programming  problem 
corresponding  to  the  above  is 

IT  T 

Problem  QP:  min  ^  x  Dx  +  c  x 

subject  to  Gx  =  h  . 

Now  we  can  obtain  very  accurately  the  following  Cholesky  decomposition 
of  D  : 

°  '  A  ’  b  ^ 


where  we  deliberately  use  A  again  to  represent  the  triangular  factor. 
If  D  is  semi-definite,  a  symmetric  permutation  of  rows  and  columns 
will  generally  be  required.  If  D  is  actually  positive  definite  then 
A  will  be  a  non-singular  triangular  matrix. 

With  the  above  notation,  it  can  be  shewn  that  the  solutions  of  both 
problems  satisfy  the  system 


f  ^ 

'  Z| 

h 

(11.1) 

I  A 

l  „T  „T  i 

r 

i 

= 

b 

V  G  A  J 

Xi 

c 

where 

c  =  0  ,  r  =  b-  Ax  for  Problem  LS, 
b  *  0  ,  r  *•  -Ax  for  Problem  QP, 


and  z  is  the  vector  of  Lagrange  multipliers.  In  [  2  ],  [3]  methods 
for  solving  such  systems  have  been  studied  in  depth.  The  method  we 
give  here  is  similar  but  more  suited  to  our  purpose.  This  method  has 
been  worked  on  independently  by  Leringe  and  Wedin  [17].  The  solution 


of  ui.l)  is  not  unique  if  the  quantity  rank^^^  is  less  than  n  , 
but  in  such  cases  we  shall  be  content  with  obtaining  one  solution  rather 
than  many.  The  important  steps  follow. 

T 

(1)  Let  be  the  orthogonal  matrix  which  reduces  G  to  triangular 

T 

form,  and  let  also  be  applied  to  A  ,  thus: 


(11.2)  ft1(GT  |  AT) 


f  R1  I  S  'N 
V  0  I  1 J 
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As  explained  earlier,  can  be  constructed  as  a  sequence  of 

Householder  t ran sformat ions,  and  the  columns  of  G  should  be 
permuted  during  the  triangularization .  This  allows  any  redundant 
constraints  in  -Jx  -  h  to  oe  detected  and  discarded. 

rp 

(2)  Let  be  the  orthoror.al  matrix  whion  reduces  '!’*  to  triangular 

form: 


Here  we  assume  for  simplicity  that  T  is  of  full  rank,  which  is 
equivalent  to  assuming  that  (11.1)  has  a  unique  solution,  and 
again  we  suppress  permutations  from  the  notation. 

(3)  The  combined  effect  of  these  decompositions  is  now  best  regarded 
as  the  application  of  an  orthogonal  similarity  transformation  to 
system  (11. l),  since  the  latter  is  clearly  equivalent  to 


The  resulting  system  consists  of  various  triangular  sub-systems 
involving  ,  Bg  ,  S  ,  and  can  easily  be  solved. 

(4)  If  desired,  the  solution  thus  obtained  can  be  improved  upon  via 
the  method  of  iterative  refinement  [ 9  ],  since  this  just  involves 
the  solution  of  system  (ll.l)  with  different  right-hand  sides,  and 
the  necessary  decompositions  are  already  available. 
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The  algorithm  Just  described  has  been  tested  on  extremely  ill-conditioned 
systems  Involving  inverse  Hilbert  matrices  of  high  order  and  with  iterative 
refinement  has  given  solutions  which  are  accurate  to  full  machine  precision. 


3^ 


12 


Positive  Definite  Programninr 


With  the  algorithm  oi'  the  previous  section  available,  we  are  now 
prepared  to  attack  the  following  more  general  programming  problems: 

Problem  LS:  min  -ib  -  Axllg 

subject  to  'J^x  ■=  , 

u,,x  >  h0  • 


IT  T 

Problem  QP:  rr.in  —  x  Dx  +  c  x 

subject  to  the  same  contraints. 


Let  G^Gg  be  of  orders  p1xn  ,  p2  x  n 
that  D  has  the  Cholesky  decomposition 
consider  problems  for  which  rank 


respectively,  and  again  suppose 
T 

A  A  .  In  this  section  we 
=  n  (which  is  most  likely 


to  be  true  with  least  squares  problems,  though  less  likely  in  QP  ) . 

In  such  cases  the  quadratic  form  is  essentially  invertible  (but  we 
emphasize  that  its  inverse  is  not  computed)  and  so  x  can  be  eliminated 
from  the  problem .  With  the  notation  of  the  preceding  section  the  steps 
are  as  follows: 


(1)  Solve  (11.1)  with  G-^h^  to  get  the  solution  x  =  xQ  ,  then  compute 
the  vector  q  =  GgXQ  -  hg  . 

(2)  If  q  >  0  then  x^  is  the  solution. 

Otherwise,  transform  the  inequality  matrix  using  from  step  (1), 
so  that 
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-T  T 

as  before  and  if  M  =  Rg  V  it  can  be  shown  that 

the  active  constraints  are  determined  by  the  following  linear 
complementarity  problem  (LCP): 

w  n  q  +  M^Mz 

(12.1) 

T 

W,  z  >  0  ,  z  w  =  0  . 

v, z  are  respectively  the  slack  variables  and  Lagrange  multipliers 
associated  with  the  inequality  constraints. 

(1+)  The  active  constraints  (for  which  =  0  in  the  solution  of 

the  LCP)  are  now  added  to  the  equalities  G^x  =  h^  and  the  final 
solution  is  obtained  from  (11.1). 

We  wish  to  focus  attention  on  the  method  by  which  the  LCP  (12.1)  is 

solved.  Cottle  and  Dantzig’s  principal  pivoting  method  (  5  ]  could  be 

T 

applied  in  a  straightforward  manner  if  MM  were  computed  explicitly, 

T 

but  for  numerical  reasons  and  because  M  M  (p2  x  p2)  could  be  very 
large,  we  avoid  this.  Rather  we  take  advantage  of  the  fact  that  no  more 
than  n  -p^  inequalities  can  be  active  at  any  one  time  and  work  with  a 
basis  M^  made  up  of  k  columns  of  M  ,  where  1  <  k  <  n  -  p^  .  The  QR 
decomposition 

«"i  -(0) 

is  maintained  for  each  basis  as  columns  of  M  are  added  to  or  deleted 
from  M^  and  as  we  knew,  Q  and  R  can  be  updated  very  quickly  each 


(3)  If 


Ci2T 


i") 
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Change.  Then  just  as  in  the  IAJ  method  i’or  linear  programming,  the  new 
basic  solution  is  obtained  not  by  updating  a  simplex  tableau  but  simply 
by  solving  the  appropriate  system  o  f  equations  using  the  available 
decomposition. 

As  an  example  ve  show  how  complementary  basic  solutions  may  be 
oaxained.  Let  the  basis  contain  k  columns  of  K  and  let 
be  the  remaining  (non-basic)  columns.  The  system  to  be  solved  is 


with  obvious  notation.  If  we  define  y  =  -M^z^  this  is  best  written  as 


and  the  solution  of  (12.2)  is  readily  obtained  from 


ZB  -B‘lu  » 


4oY 

\  Of  }  n-px-k 


The  blocking  variable  when  a  non-basic  variable  is  increased  can  be 
found  from  the  solution  of  the  same  set  of  equations  with  the  appropriate 
right-hand  side.  It  is  worth  noting  that  the  equations  can  be  simplified 
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if  the  basis  is  square  (i.e.,  if  there  are  as  many  constraints  active 
as  there  are  free  variables) .  Since  it  seems  very  common  for  the  basis 
to  fill  up  during  the  iterations  (even  if  the  final  solution  does  not 
have  a  full  set  of  constraints)  it  is  worth  treating  a  fiill  basis 
specially. 

Almost-canplementary  solutions  can  be  obtained  in  similar  fashion 
(with  somewhat  more  work  required  as  the  system  is  then  not  quite  so 
symmetric) .  Thus  an  algorithm  such  aB  Cottle  and  Dantzig'a  can  be 
implemented  using  these  techniques,  and  convergence  is  thereby  guaranteed. 

Of  special  interest,  however,  is  the  following  unpublished  and 
apparently  novel  idea  due  to  Yonathan  Bard,  with  vhoBe  permission  we 
report  the  results  he  has  obtained.  Almost -complementary  bases  are 
never  allowed  to  occur;  instead,  if  a  basic  variable  is  negative, 
then  it  is  replaced  by  its  complement  regardless  of  the  effect  on  the 
other  basic  variables.  Bard  has  tried  this  method  (carried  to  convergence) 
on  hundreds  of  problems  of  the  form  v  =  q+Mz  and  cycling  has  never 
occurred  when  the  most  negative  element  of  q  is  chosen.  In  a  series 
of  tests  on  100  random  matrices  of  orders  between  2  and  20  , 
principal  pivoting  required  a  total  of  537  pivots  whereas  the 
Cottle-Dantzig  algorithm  required  689  . 

The  present  authors1  experience  with  fewer  but  larger  problems 
confirms  the  above  observation  that  convergence  does  actually  occur  and 
usually  after  a  small  number  of  iterations.  Since  the  idea  eliminates 
all  work  other  than  ccmputstion  of  complement  butt  solutions  it  is 
particularly  suited  to  the  techniques  of  this  section.  At  worst  it  should 
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be  used  as  a  starting  procedure  to  find  u  close -to-optimal  basis  quickly, 
and  at  best  If  the  conjecture  can  be  proven  that  it  will  always  converge, 
then  a  lot  of  computer  time  could  be  saved  in  the  future. 

[It  has  since  been  learned  that  Bard  applied  the  principal-pivoting 
rule  to  LCP's  of  the  somewhat  special  form  in  which 

M  *  P*P,  q  «  Pxp 

for  some  P,  p.  Problems  of  this  form  have  been  studied  by  Zoutendijk 
In  [31,321  where  several  pivot  selection  rules  are  discussed.  Finite- 
ness  is  proven  for  one  rule,  but  simpler  methods  (such  as  Bard's)  are 
recommended  in  practice  for  efficiency. 

The  question  of  finiteness  for  the  more  general  ICP  remains  open, 
and  it  is  likely  that  somewhat  more  sophisticated  rules  (e.g. ,  Cottle 
and  Dantzig)  will  be  required. 
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We  nov  consider  the  more  general  problem  in  which  the  rejik  of  the 
quadratic  form  combined  with  the  equality  constraints  may  be  less  than 
The  method  we  propose  is  conceptually  as  simple  as  it  is  stable.  It  is 
analogous  to  the  revised  simplex  method  for  linear  programming  in  that 
the  essential  steps  to  be  implemented  are  as  follows: 

(1)  Find  the  current  basic  solution  from  a  certain  system  of  equations 
for  which  a  decomposition  is  available. 

(2)  Determine  according  to  a  certain  set  of  rules  What  modifications 
should  be  made  to  the  system  to  obtain  a  new  basis. 

(3)  If  necessary,  update  the  decomposition  and  return  to  step  (1). 

Thus,  suppose  that  the  current  basis  contains  G^x  =  hg  as  active 
constraints.  As  in  (11.1)  the  corresponding  basic  solution  is  then 
given  by 

f  °B 

(13.1)  I  A 

U  i 

and 

(!3.2)  VB  =  *  5BX  ' 

(Here,  Gflx  >  ^  are  the  currently  inactive  constraints,  wg  the 
corresponding  slack  variables,  and  the  Lagrange  multipliers  or  dual 
variables  associated  with  the  active  constraints.)  The  elements  of  z 

B 


corresponding  to  any  equality  constraints  may  be  either  positive  or 
negative  and  need  never  be  looked  at.  Ignoring  tnese,  the  basic  solution 
above  is  optimal  if  and  only  if 

>  0  and  w_,  >  0 

B  -  3  — 


A  "QP  algorithm”  is  r.ow  to  be  regarded  as  the  "certain  set  of  rules" 


mentioned  in  step  {2}  whereby  z and  possibly  other  information  are 
used  to  determine  which  constraints  should  be  added  to  or  dropped  from 
The  efficiency  of  the  method  will  depend  on  the  speed  with,  which  this 
decision  can  be  made  and  on  the  efficiency  with  which  the  decomposition 
of  (13.1)  can  be  updated. 

Once  again  the  most  promising  pivot- select ion  rule  is  that  of  Bard, 
as  discussed  in  the  previous  section.  The  general  idea  in  this  context 
is  as  follows: 


(a)  Find  •  =  min  w.  ,  z  =  min  z.  from  those  eligible 

1  p  1 

elements  of  wD,z„  . 

d  a 

(b)  If  va  <  0  ,  constraint  a  could  be  added. 

(c)  If  Zp<0,  constraint  p  could  be  dropped. 

(u;  If  there  are  already  n  constraints  active  and  w  <  0  , 

OL 

corstraint  Q  could  replace  constraint  0  . 


We  do  not  consider  here  the  question  of  convergence,  but  as  already  stated, 
this  type  of  rule  has  been  found  to  work. 

The  problem  of  updating  the  requisite  decompositions  is  more  relevant 
at  present.  We  discuss  this  and  other  points  briefly. 


(1)  The  matrices  Q  ,1^  p  Equation  (11.2)  can  be  updated  efficiently 
using  the  methods  of  'ection  8. 

(2)  Q2>R2  Gained  from  the  matrix  T  in  Equation  (11.3)  unfortunately 
cannot  he  updated,  but  the  work  needed  to  recompute  then  might  often 
be  very  small,  for  the  following  reasons: 

(a)  In  Problem  LS,  a  preliminary  triangularization  of  A  (tutu) 

can  be  applied  to  obtain  an  equivalent  problem  for  which  m  <  n  . 
The  Choiesky  factor  of  D  in  Problem  QP  already  has  this  property. 

(b)  If  there  are  many  constraints  active  (up  to  n)  then  T  has 
very  few  rows. 

(c)  If  the  rank  of  the  system  is  low  (relative  to  n)  then  T 
has  very  few  columns. 

(3)  Hence  the  method  is  very  efficient  if  close  to  n  constraints  are 
active  each  iteration,  as  should  often  be  the  case.  It  also  has  the 
property,  dong  with  Beale's  algorithm  [l ],  of  being  most  efficient 
for  problems  of  low  rank. 

(b)  The  procedure  can  be  initiated  with  any  specified  set  of  constraints 
in  the  first  basis,  and  an  initial  estimate  of  x  is  not  required. 

(5)  Any  number  of  constraints  can  be  handled,  in  the  same  way  that  the 
revised  simplex  method  can  deal  with  any  number  of  variables. 

(6)  If  D  =  0  the  problem  is  a  linear  program  and  only  bases  containing 
n  constraints  need  be  considered.  The  method  reduces  to  something 
like  a  self- dual  simplex  algorithm. 


k2 


Finally  we  note  that  with  ceKi-loflni.lt:  problems  it  is  possible 
for  some  basic  system  (1J.1)  to  be  a  Smaller.  if  tr.ere  arc  any  solutions 
at  all  then  there  are  many  (tills  will  always  be  the  case  with  low  rank 
least  squares  problems)  but  this  does  not.  matter,  since  :•  is  still 
uniquely  determined.  However,  a  .low  rank  quaarar.ic  program  might  be 
unbounded,  ar.d  this  is  manifested  by  a  singular  system  (1(1 .  l)  proving 
to  be  inconsistent.  In  general,  inis  just  means  tnat  there  arc  not  yet 
enough  constraints  in  the  basis,  so  that  trouble  can  usually  be  avoided 
by  .initialising  the  procedure  wit::  a  full  set,  of  constraints. 
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1^.  The  Singular  Value  Decomposition 

Let  A  be  a  real,  sn  matrix  (for  rotational  convenience  we 
assume  that  in  >  n)  .  It  Is  well  known  (cf.  [  ])  that 

'.n.l)  A  =  :  z/ 

wnere  U,  V  arc-  orthogonal  matrices  and 


'■  consists  of  the  orthoncrmaiized  eigenvectors  of  AAT,  and 

T 

V  consists  of  the  orthonormalized  eigenvectors  of  A  A  .  The 

diagonal  elements  of  Z  are  the  non-negative  square  roots  of  the 
T 

eigenvalues  of  A  A  ;  they  are  called  singular  values  or  principal  values 
of  A  .  We  assume 

°1  —  —  *  *  *  —  5n  —  0  * 

Thus  if  rank(A)  =  r  ,  =  cp+2  =  ”*  =  en  =  0  .  The  decomposition 

( l1*  •  1)  is  called  the  singular  value  decomposition  (SVD). 
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An  n  x  m  matrix  X  is  said  to  be  the  pseudo -inverse  of  an  m  x  n 


matrix  A  if  X  satisfies  the  following  four  properties: 

(i)  AXA  =  A  ,  (ii)  XAX  ■=  X  ,  (iii)  (XA)T  =  XA  ,  (iv)  (AX)T  =  /iX 

V.'e  denote  the  pseudo-inverse  by  A+  .  It  can  be  shown  that  A+  can 
always  be  determined  and  is  unique  (cf.  L  2li) •  It  is  easy  to  verify 

+  T 

that  A  =  VAU  where  A  is  the  nxm  matrix 

A  -  diagta^a^1, .  .  ..,0]  .  There  are  many  applications  of 

the  S7B  in  least  squares  problems  (cf.  [11]). 

The  SVD  of  an  arbitrary  matrix  is  calculated  in  the  following  way. 
First,  a  sequence  of  Householder  transformations 
is  constructed  so  that 


PnPn-l‘*‘PlAQiq2' 


’Vi  s  AQ  -  J 


and  J  is  an  m  x  n  bi -diagonal  matrix  of  the  form 


The  singular  values  of  J  are  the  same  as  those  of  A  . 

Next  the  SVD  of  J  is  computed  by  an  algorithm  given  in  [ll  ] .  The 
algorithm  is  based  on  the  highly  effective  Q.R  algorithm  of  Francis  [10]  for 
computing  eigenvalues.  If  the  SVD  of  J  -  X£YT  then  A  =  PX2YTQT  so 
that  U  =  HC  ,  V  *  QY  . 
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15-  Nonlinear  Least  Squares 

Consider  the  nonlinear  transformation  F(x)  -  y  where  . 
and  ycEm  with  n  <  m  .  We  wish  to  consider  the  following  problem: 

min  !'b-F(x)!|2 

subject  to 

(15-1)  Gx  =  h  f 


where  0  is  a  p  x  n  matrix  of  rank  p  and  heE  .  A  very  effective 

P 

algorithm  for  solving  such  problems  is  a  variant  of  the  Levenberg-Marquardt 
algorithm  [18, 19];  in  this  section  we  consider  some  of  the  details  of  the 
numerical  calculation.  Further  extensions  of  the  algorithm  are  given 
by  Shanno  [25]  and  Meyer  [20]. 

Let  us  assume  that  we  have  tun  approximation  which  satisfies 

the  relation  Gx^  »  h  .  Then  at  each  stage  of  the  iteration  we 
determine  so  that 


(15-2) 

and 

(15.3) 


x<k+D  =  x<k>  +  6^ 

G6^  =  0  . 


Again  as  in  Section  11,  we  write  Q  ^  k J  where  is  the  product 
of  p  Householder  transformations  and  R  is  an  upper  triangular  matrix. 
Let 


(15.k) 


Then  from  (15.3),  we  see  that  _  0  . 

For  notational  convenience,  let  us  drop  the  superscript  k  ; 
we  write  as  x^  ar.d  as  x  . 

In  the  Levenberg-Marquardt  algor  itlim  one  Jete  rainy  a  the  vector  £, 
so  that 

(15-5)  Mr  -  J^|S^  +  xlifif  =  min. 

where 


r  =  b  -  F(xq)  , 


J  is  the  Jacobian  evaluated  at  x  ,  and  \  is  an  arbitrary  non-negative 

U 

parameter.  From  (15A)j  we  see  that  (15-5)  is  equivalent  to  determining  T] 
so  that 


(15.6) 


{ 


subject  to 


+ 

I 


X(||S|la  +  111  III  )= 
=  0  . 


T 

Now  let  us  write  =  [M,N]  where  N  consists  of  the  last  n-p 

T 

columns  of  JQ^  .  Then  (15*6)  is  equivalent  to  finding  1]  so  that 

*0l)  =  ||r-  Wlllg  +  Xlhllg  =  min. 

Consider  the  SVD  of  N  ;  namely 
N  -  UIV2  . 


U7 


Then 


(i5.7)  *(n) 


-  !'e  dil +  >•!!  5II2 

where 

s  =  UTr  ,  £  =  VTT1  . 

Writing  out  (15-7)  explicitly,  we  have 

,(t)  ‘  &(*i  -  °jt/  * 1  %  <«/ 

where  p  is  the  rank  of  N  .  (Note  p  nay  change  from  iteration  to 
iteration.)  Then 

A  • 

$(5)  -  min 

when 


;  .  _f£L 

X  +  cr. 


for  j  =  1,2, ..  ,,p  , 


=  0 


for  j  >  p 


ar.d  hence 


where  v  is  the  j-th  column  of  V  .  Thus 
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we  assume  s.  f  o  for  j  -  1,2,  ...,p  .  By  repeated  use  of  the 
relationship 


det(z  w)  *  det(X)  det(W-ZX_;1Y)  if  det(X)  /4  0 

we  can  show  that  (15*9)  is  equivalent  to 
(15-10)  det((n  +  II)2  -  uuT)  =  0 

which  has  2p  roots;  it  can  be  shown  that  we  need  the  largest  real 
root,  which  we  denote  by  X  ([8])-  Let 


2  2  2  2  2 
and  assume  that  >  o2  >  ...  >  a  >  0  .  Note  r(0)  =  g  -  a  >  0  , 

2 

and  r(x)  -*  -a  as  \  ®  ,  so  that  0  <  X  <  03  ®nd  it  is  the  only 

* 

root  in  that  interval.  We  seek  a  more  precise  upper  bound  for  X  • 


50 


From  (1^.10)  we  see,  using  a  Rayleigh  quotient  argument, that 
\  <  max  ( -y*  Cl  y  * 

"•  II  »J 

'iy  l2=x 

A  short  manipulation  then  shows  that 


Thus  ,  wo  wi$h  to  find  a  root  of  (15-10)  which  lies  in  the  interval 
hiven  by  (I5.U).  Hole  that  the  deterainantal  equation  (I5.I0) 
involves  a  diagonal  matrix  plus  a  matrix  of  rank  one.  In  the  next 
section  we  shall  describe  ar,  algor  it  San  for  solving  such  problems. 


16.  Modified  Elgensystans 

As  was  pointed  out  in  Section  15,  it  is  sometimes  desirable  to 

determine  some  eigenvalues  of  a  diagonal  matrix  which  is  modified  by 

a  matrix  of  rank  one.  Also,  Powell  [23]  has  recently  proposed  a 

minimization  algorithm  which  requires  the  eigensystem  of  a  matrix  after 

a  rank  one  modification.  In  this  section,  we  give  an  algorithm  for 

determining  in  0(n  )  numerical  operations  some  or  all  of  the  eigenvalues 

T 

and  eigenvectors  of  D  +  ouu  where  D  =  diag(di)  is  a  diagonal  matrix 

of  order  n  and  ueE  . 

n 

T 

Let  C  =  D+tfuu  ;  we  denote  the  eigenvaluesof  C  by  . .  ,,\n 

and  we  assume  >  xi+1  and  >  di+1  .  It  can  be  shown  (cf.  [30 J) 
that 

(1)  If  o  >  0  ,  dj^  +  CT^u  >  >  dx  ,  di_1  >  Xi  >  d±  (i  =  2,  ...,n)  , 

(2)  If  o<0,  di>^i>di_i  (i  =  1,2, . .  .,n-l)  ,  dn  >  3^  >  dn  +  ffuTu  . 

Thus  we  have  precise  bounds  on  each  of  the  eigenvalues  of  the  modified 
matrix. 

Let  K  be  a  bi-diagonal  matrix  of  the  form 


and  let  M  =  diag(^i)  .  Then 


is  a  symmetric,  tri-diagonal  matrix. 

Consider  the  matrix  equation 

(16.2)  (D  +  ctuu  )x  =  ?^X 

Multiplying  (16.2)  on  the  left  by  K  ,  ve  have 

K(D  +  ouuT)  kVTx  =  \  K  KTK'Tx 
or 

(16.3)  (KDif +  aKuuTKT)y  =  UCK^y 

where  x  =  K?y  .  Let  us  assume  that  we  have  re-ordered  the  elements  of  u 
so  that 

\  -  Ug  -  ...  =  =  0  and  0  <  |up|  <  <  ...  <  |uj  . 
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Now  it  is  possible  to  determine  the  elements  of  K  so  that 

(USA)  Ku  = 

Specifically, 

=  0  (i  =  1,2,  ...,p-l)  , 

ri  =  "ui/ui+l  (i  =  ••*>«)  , 

and  we  note  that  |i\  |  <1  .  (This  device  of  using  a  bi-diagonal  matrix 
for  annihilating  n-1  elements  of  a  vector  has  been  used  by  Bj&rck 
and  Pereyra  ( It  ]  for  inverting  Vandermonde  matrices.)  Therefore,  if  Ku 
satisfies  (l6.lt),  we  see  from  (l6.1)  that  KDJ?  +  crKuu^K1'  is  a 
tri-diagonal  matrix  and  similarly  KK1  is  a  tri-diagonal  matrix.  Thus 
we  have  a  problem  of  the  form 

Ay  =  \By 

where  A  and  B  are  symmetric,  tri-diagonal  matrices  and  B  is  positive 
definite. 

In  [22],  Peters  and  Wilkinson  show  how  linear  interpolation  may 
be  used  effectively  for  computing  the  eigenvalues  for  such  matrices 
when  the  eigenvalues  are  isolated.  The  algorithm  makes  use  of  the  value 
of  det(A-XB)  .  When  A  and  B  are  tri-diagonal,  it  is  very  simple 
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to  evaluate  det(A  -  xB)  for  arbitrary'  \  .  Once  the  eigenvalues  are 
computed  it  is  easy  to  compute  the  eigenvectors  by  inverse  iteration. 

In  Section  15>  we  showed  it  was  necessary  to  compute  a  parameter 

* 

X  which  satisfied  the  equation 

(16.5)  dert((n  +  XI)2  -  uu?)  =  0  . 

Again  we  can  determine  K  so  that  Ku  satisfies  (16.4)  and  hence  (16. 5) 
is  equivalent  to 

(16.6)  det (K(Q  +  XI)2KT  -  KuuTl?)  *=  0  . 

The  matrix  G(x)  =  K(Q  +  Xl)^KT  -  Ku^K1  is  tri-diagonal  so  that  it  is 

easy  to  evaluate  G(x)  and  det  3(x)  .  Since  we  have  an  upper  and 

lower  bound  on  X  >  it  is  possible  to  use  linear  interpolation  to 
¥• 

find  X  »  even  though  G(x)  is  quadratic 'in  \  .  Numerical  experiments 

have  indicated  it  is  best  to  compute  G(X)  -  K(G  +  Xl)2KT  -  Kuu^K1 

* 

for  each  approximate  value  of  \  rather  than  computing 
G(X)  =  (KQ2  K1  -  KuuV)  +  GXKOK1  +  \2  Kt?  . 

The  device  of  changing  modified  eigensystems  to  tri-diagonal 
matrices  and  then  using  lir.ear  interpolation  for  finding  the  roots  can 
be  extended  to  matrices  of  the  form 


Again  we  choose  K  so  that  Ku  satisfies  (16.4)  and  thus  obtain  the 
eigenvalue  problem  Ay  =  XBy  where 
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I 

Ht1  0  \ 

°  I1/ 

so  that  A  and  B  are  both  tri -diagonal  and  B  is  positive  definite. 
Bounds  far  the  eigenvalues  of  C  can  easily  be  established  in  terms  of 
the  eigenvalues  of  D  and  hence  the  linear  interpolation  algorithm 
may  be  used  for  determining  the  eigenvalues  of  C  . 


I 

I 

I 

I 
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